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ABSTRACT 

2D axis-symmetric hydrodynamical simulations are presented which explore the inter- 
action of stellar and disk winds with surrounding infalling cloud material. The star, 
and its accompanying disk, blow winds inside a cavity cleared out by an earlier jet. The 
collision of the winds with their surroundings generates shock heated plasma which 
reaches temperatures up to ~ 10^ K. Attenuated X-ray spectra are calculated from 
solving the equation of radiative transfer along lines-of-sight. This process is repeated 
at various epochs throughout the simulations to examine the evolution of the intrinsic 
and attenuated flux. We find that the dynamic nature of the wind-cavity interaction 
fuels intrinsic variability in the observed emission on timescales of several hundred 
years. This is principally due to variations in the position of the reverse shock which 
is infiuenced by changes in the shape of the cavity wall. The collision of the winds 
with the cavity wall can cause clumps of cloud material to be stripped away. Mixing 
of these clumps into the winds mass-loads the fiow and enhances the X-ray emission 
measure. The position and shape of the reverse shock plays a key role in determining 
the strength and hardness of the X-ray emission. In some models the reverse shock is 
oblique to much of the stellar and disk outfiows, whereas in others it is closely normal 
over a wide range of polar angles. For reasonable stellar and disk wind parameters the 
integrated count rate and spatial extent of the intensity peak for X-ray emission agree 
with Chandra observations of the deeply embedded MYSOs S106 IRS4, Mon R2 IRS3 
A, and AFGL 2591. 

The evolution of the cavity is heavily dependent on the ratio of the infiow and 
outfiow ram pressures. The cavity closes up if the infiow is too strong, and rapidly 
widens if the outfiowing winds are too strong. The velocity shear between the respec- 
tive fiows creates Kelvin-Helmholtz (KH) instabilities which corrugate the surface of 
the cavity. Rayleigh- Taylor- like instabilities also occur when the cavity wall is pushed 
forcefully backwards by strong outfiows. The opening angle of the cavity plays a sig- 
nificant role and we find that for collimation factors in agreement with those observed 
for bi-polar jets around massive young stellar objects (MYSOs), a reverse shock is 
estabhshed within < 500 au of the star. 

Key words: X-rays:stars - stars:winds, outfiows - stars: formation - stars:early-type 
- hydrodynamics 



1 INTRODUCTION 

Observational and theoretical advances have provided in- 
creasing evidence that massive star formation is not merely a 
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scaled up version of low mass star formation. The crucial dif- 
ference arises from the promine nt role of r adiation pressure 
in massive star formation (^Zinnecker York e 2007, and ref- 
erences there-in). Theoretical models find that the radiation 
pressure from the core is too great for s pherical accretion to 
form stars above a mass of - 10 Mo (IWolfire Cassinellil 
ll987l V This problem di sappears if one considers disk ac - 
cretion (lYorke Sonnha lter 20^; iKrumholz et al.1 l2009l ). 
Observations of young massive stars support this scenario 
as they show the presence of a dense equatorial disk and 
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on^oin^ accretion (|Patel et all l2005l : iBeltran et alJ l2006l : 
iTorrelles et"aIll2QQ7l ). There are, however, some similarities 
in the formation processes of both lo w and hi^h m ass stars. 
For instance, both invo l ve outflows (|Garav Lizano 19991 : 
iReipurth Ballvl l2Q0ll : iBaneriee Pudritd l2QQ6l . l2007f ). 
Bi-polar cavities a round YSOs are com monly observed 
in star formation (jGarav Lizanol Il999l). with collima- 
tion factors of ^ 2-10 for MYSOs (iBeuther et alJ l20Q2l : 
IPavis et al.ll20o3 V The cavities may be more strongly col- 
limated at large distances (due to their formation by jets), 
but less so at the base due to the action of stellar winds 
(jShepherdl [2005). For low -mass stars the observed cavity 
morphology can be qualitatively reproduced by the injec- 
tion of an outflow into a non-spherical density distribution 
(iDelamarter et al.1 |200Q|: iLee et alJ l200ll: Iwilkin Stahlerl 



l2QQ3l : iMendoza et al.ll2004l : [Cunningham et al. "2005V Mag- 
netic field s increase the collimation (Gardiner et al. l2003l : 
IShang et alJ l20Q6l ). Numerical simulations of high-mass 
star formation show that outflows are a by-product of 
the collapse process which forms the massive star, and 
may be due to radiation press ure and m agnetic fields 
(lYorke Sonnha lter 2002: .Baneriee Pudr itz 2007). 

The widths of IR recombination line emission observed 
from MYSOs indicates the presence of dense o utflows with 
veloc i ties ranging from 100 to > 340 kms~^ (|Drew et al.l 
Il993l : iBunn et all Il995l l. Further confirmation of outflows 
has come from high angular resol ution radio observations 
(e.g. iHoare et"^ ] |l994l : lHoarj|2006l ). One explanation would 
be a disk wind generated when UV flux from the star is 
absorbed and re-emitted by the material at the surface of 
the disk. 

Radiation-driven accretion disk winds ha ve been pre- 
viousl y exp lored using h ydrod ynamic models. IProga et al.l 
(|l998l ) and IPrew et al.1 (jl998l ) found that for a star-disk 
system where the 10 M© star and it's accompanying disk 
were equally luminous (i.e. L*/Ldisk = 1), a strong steady 
disk wind could be produced. They found that there was a 
transition from a fast, low density polar outflow to a higher 
density outflow at a polar angle of ^ ^ 40 — 60° . The stellar 
wind had a velocity of 2000 kms~^, and the disk wind ve- 
locity ranged from 400 km s~^ at high latitudes to tens of 
kms~^ at low latitudes. They also found a density contrast 
of 10^ — 10^ between the stellar and disk winds, the latter 
being much denser. The opening angles of the stellar and 
disk winds are dependent on a number of parameters in the 
disk wind model: the luminosity ratio between the star and 
the disk, the forc e parameter used to calculate the radiative 
driving (jProga et al. 199^), and the magnetic pressure in 
the wind (jProgal 120031 ). However, the assumption of strong 
magnetic fields in MYSO disks is speculative. The equatorial 
flow detected around S106 I RS4 (or IR) show s a mass- loss 
rate of M ^ lO~^M0yr"^ (IPrew et al.lll993l : iHoare et al.1 
Il994l ). a factor of ^ 100 higher than the simulations. Such a 
difference c ould be accounted for by clumping in the wind. 
ISim et al.1 (l2005l ) used a radiative transfer model inspired 
by disk wind simulations to model the HI line emission 
from MYSOs. They adopted a disk wind which was largely 
equatorial and had a velocity which declined from ^ 850 
to 300kms~^ in the angular range 82 — 89°. This model 
proved to be successful in reproducing qualitative features 
of Brackett and Pfund series line emission observations. 

X-rays have been detected from deeply embedded 



MYSOs in star forming regions by the Chandra X- 
ray observatory (hereaf ter Chandra) (e.o:. iBroos et al.ll2007l : 
IWang et al] 1200/1 . l2009l ). An initial observation of Mon R2 
detected X-ray emission from the intermediate mass ob- 
jects consistent with the hard and highly time variable X- 
ray emission caus ed by magnetic flaring activity between 
the star/disk (Ko hno et ahlfjoO^ . Further analysis of the 
same data, coupled with high resolution near-IR interfer- 
ometry identified the s eparate constituents of Mon R2 IRS3. 
IPreibisch et al.l (|2002h found that the X-ray emission from 
IRS3 A and C, with a count rate of 0.166 ± 0.041 ks"^ 
for the former, could not be explained by the standard 
scenario for massive stars (i.e. wind embedded shocks pro- 
du ced by instabilities i nherent in radiatively-driven winds - 
see lOwocki et al1ll988l ), yet the estimated stellar masses of 
these objects implies they will have radiative outer envelopes 
which poses problems for the generation of X-rays through 
magnetic sta r /disk interactions. In their observations of the 
S106 region, iGiardino et al.l (|2004l ) suggested that it may 
not be necessary to employ magnetic fields to explain the 
X-ray emission characteristics of the massive central object 
S106 IRS4 where the count rate was 0.30 =b 0.11 ks"^ 

A potential source of X-rays which has not been con- 
sidered before is the collision between the stellar and disk 
winds and the infalling envelope. To explore this interaction 
we perform the first detailed modelling of an MYSO embed- 
ded at the centre of a bi-polar cavity. The star and disk drive 
supersonic outflows, and we examine how these outflows in- 
fluence the evolution of the infalling envelope. We find that 
the wind-cavity interaction produces a reverse shock where 
gas is heated to X-ray emitting temperatures. Our models 
are able to reproduce the observed count rates, although 
there is a strong sensitivity to the inclination angle of the 
observer. 

This paper is organised as follows: § [2] gives a descrip- 
tion of the model used in the simulations. The results of the 
simulations and the observable X-ray emission are presented 
in §[3l In §[4lwe report on a recent observation of the MYSO 
AFGL 2591 with Chandra then perform fits to this observa- 
tion and those of two other candidate objects, § [5] includes 
discussion, and in § [6] we summarize the conclusions from 
the current work and suggest possible future directions. 



2 THE MODEL 

Our model is of a recently formed massive star with an ac- 
companying accretion disk. The star and disk are situated at 
the centre of a flattened envelope with a pre-existing bi-polar 
cavity evacuated by earlier jet activity. The morphology of 
the cavity is prescribed by a simple analytic equation. A 
schematic of the simulation domain is shown in Fig. [1] 

The cavity intersects the disk plane at the cavity radius, 
Tcav In essence, the presence of a disk is implied through 
the presence of the disk wind, yet we do not include any 
parameters for the disk itself. We simply note that the disk 
wind is la unched from the inner regions of the disk (i.e. 
r < 10 R.. lDrew et~alll 19981 ). The star and disk are situated 
in the first cell within the simulation domain. Initially the 
cavity is filled with unshocked stellar and disk wind material, 
with the density and velocity a function of the polar angle, 0. 
The material outside of the cavity is infalling, and represents 
molecular gas which is still accreting onto the disk. 
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Figure 1. Schematic of the simulation domain. The left-hand 
panel shows a zoom in of the region close to the star and disk. 
The right panel shows the entire simulation domain, with dotted- 
line representing the edge of the grid. The cross-hatched quadrant 
in the left-hand panel represents the region where the winds are 
mapped onto the grid and the adjacent rectangular region repre- 
sents the additional cells used to ensure the angular dependence 
of the wind is resolved (see ^ 12. 2|) . 



Features of the model, details of the hydrodynamics 
code used, and a description of the X-ray calculations are 
now given. 



2.1 The infalling cloud 

For the cloud material we use the solutions for the collapse of 
a slowly rotating isotherm al sphere with a single point source 
of gravity at the centre (lUlrichI [19761: ICassen MoosmanI 
Il98ll : IChevalierl [19831 : iTerebev et al.lll984l ). With conserva- 
tion of angular momentum for infalling material the tra- 
jectories are parabolic and are described by the streamline 
equation (in polar coordinates). 



^ R (1 



Mo)Mo 



(1) 



where /x = cos^, /io = cos^o, and R — + Tc is the 
centrifugal radius of the cloud and Oo is the initial polar 
angle. The components of the infall velocity are given by. 
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The density for the cloud material is then found by integrat- 
ing the mass flux along streamlines. 



[1 + (3m§ - 1)C]" 



(3) 



The infalling cloud material is assumed to be molecular 
gas at a temperature of 100 K. 



2.2 The winds 

Following the dete c tion of equa t orial w inds from MYSOs 
(iDrew et al l Il993l : iHoare et all Il994l : iBunn et al.1 Il995l : 
lHoarell2006l i we choose to incorporate a disk wind into our 
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Figure 2. Variation of the angle dependent mass-loss rate, Mq 
(upper panel) and wind velocity, vq (lower panel) as a function 
of polar angle, ^, for model Rl. 



model. The wind geometry we adopt is aimed at represent- 
ing the characteristics of disk wind models. In particular 
w e base our w i nd ge ometry on the hydrodynamical model 
of iDrew et al.l (|l998ll and th e angle dependent model con- 
structed bv lSim et al.l (|2005l ) . In these works the disk wind 
is assumed to be launched from the inner region (< 10 i?*) 
of the accretion disk. A considerable difficulty in construct- 
ing this model is the range of scales across which impor- 
tant physics occurs. The cavity has dimensions of the or- 
der of 10^^ cm, whereas the acceleration region of the winds 
is < 10^^ cm. Therefore, to represent the winds we have 
adopted a simple latitude dependent prescription, where the 
winds are mapped on with terminal velocities, which aims to 
capture the important features of a stellar/disk wind com- 
bination. The mass-loss rate is then, 



Me 
A 



Me 
2 



[1 + A+(1- A) tanh((^-^s)x)]: 



(4) 



Me is the mass-loss rate of the outflow (stellar and disk 
winds) as a function of polar angle, is the transition angle 
between the stellar and disk winds, and x is a steepening 
factor to sharpen the transition from the stellar wind to the 
disk wind (consequently, the frictional heating between the 
two flows is linked to this parameter) . Ms and Md would be 
the mass- loss rates of the stellar and disk wind, respectively. 
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if the mass-loss was isotropic and spherical. Therefore, the 
actual net mass-loss rates of the stellar and disk winds can 
be found by integrating Me over the appropriate angular 
range. We set % = 5. The wind velocity as a function of 
polar angle, ve, is. 



i|i[l+7 + (l-27)tanh(((?-(?s)x)] 
ve = { i^[l+7 + (2-7)(^)] 




o<0<es 

Os<0<Od 

eo<0< 7T/2 



where j = vs/vu, vs and vu are the terminal velocities of the 
stellar and disk winds respectively, and is the transition 
angle above which there is no wind blown from the disk. The 
wind density as a function of polar angle is then calculated 
from. 



Pe 



Me 



ve 



(6) 



The angle dependent mass-loss rate and velocity profiles 
for model Rl are shown in Fig. [2] The winds are mapped into 
cells within a small radius of the grid origin 1 x 10^^ cm) 
provided that prior to the mapping a cell contains wind ma- 
terial (i.e. not inflowing cloud material). When necessary, 
additional cells in the first row above the disk plane (the r- 
axis in the model) are used to map on the winds (i.e. when 
On ^ 89°). This is to ensure that the angular dependence 
of the disk wind is resolved sufficiently. The stellar and disk 
winds are switched on at t = 0. The unshocked winds are 
assumed to be ionized gas at a temperature of 10^ K. 



2.3 The cavity 

In the model the initially un-shocked stellar/disk winds are 
separated from the infalling cloud material by a cavity wall. 
The shape and position of the cavity wall is described by a 
simple analyt ical prescription relating z and r on its surface 
([Alvarez et al . 2004a), 



\ I cav / 



where k ■ 



r 

acosuj 



since-) 



(7) 
(8) 



(jj is the half opening angle of the cavity, (3 controls the initial 
curvature of the cavity wall (/3 = 1 or 2 corresponds to a 
cavity wall with either a straight or parabolic shape), and 
a is a characteristic scale size of the cavity (we assume that 
a = 8 X 10^^ cm, approximately in agreement with length 
scales for bi-polar outflows around MYSOs). It should be 
noted that Eq.[8]is only suitable for a; ^ 3°. 



2.4 The hydrodynamical code 

To perform the simulations we use the hydrodynamical 
code VH-1 (|Blondin et al.1 Il990l). The code utilizes the 
Piecewise Parabolic Method ( Colella Woodward! Il984l ) 
to solve the gas dynamics equations at cell boundaries. 
The version of the code we use in this work has under- 
;one some modifjcations w hich include: radiative cooling 
Strickland Blondin|[T995l ). a direct Eulerian solver, and 
portability to multi-processor machines using the message 
passing interface (MPI) via the grid parallelization scheme 



described in lSaxton et al.l (|2005l V We include radiative cool- 
ing for optically thin gas above T = 10^ K, and assume that 
molecular gas below this temperature instantaneously cools 
to a temperature of T = 10^ K. Gravity and rotation due 
to a central point source are incorporated through effective 
potential terms. The code also includes advected scalar fluid 
variables (dyes) which are used to trace the position of the 
'^^^parate components of the flow, i.e. stellar wind, disk wind, 
or cloud material. 

We use a 2D, cylindrically symmetric, fixed grid with 
constant resolution throughout the simulations of = = 
8.3 X 10^^ cm, which for the grid of the fiducial model (Rl in 
Table [U equates to r x z = (5 x 10^^ cm) x (8 x 10^^ cm) = 
600 X 960 cells. The boundary conditions (BCs) are set as 
follows: the r = boundary uses a symmetric BC, the z = 
boundary uses an outflow BC (to allow gas to flow onto a 
hypothetical disk). The z = Zmax and r = rmax boundaries 
use either a zero- gradient outflow BC if the cell adjacent to 
the boundary is wind material or a fixed inflow condition for 
cloud material if the adjacent cell contains cloud material. 
In the latter case Eqs.[2]and[3]are used to populate the ghost 
cells. 



2.5 X-ray emission 

To allow a comparison to be made between Chandra X-ray 
observations of MYSOs and our simulations we calculate 
attenuated X-ray fluxes. Radiative transfer calculations are 
performed on the 2D m odels using the method described in 
iDoughertv et all (l2003l ). This involves calculatine: emission 
and absorption coefficients for each cell on the grid and inte- 
grating the radiative transfer equation. For the emissivities 
we use values for optically thin gas in collisional ionization 
equilibrium obtained from look-up tables calcu lated from 
the MEKAL plasma code (|Liedahl et al.l [l995l . and refer- 
ences there-in) containing 200 logarithmically spaced energy 
bins in the range 0.1-10 keV, and 101 logarithmically spaced 
bins in the range 10^ — 10^ K. Solar abundances are assumed. 
The advected fluid variables (dyes) are used to identify the 
X-ray emission contributions from the stellar and disk winds 
and the heated cloud material. 

For lines-of-sight which exit the grid while in infalling 
cloud material we calculate an additional column density be- 
yond the grid boundary by numerically integrating Eq.[3] X- 
rays which pass through the disk plane within the cavity ra- 
dius, Tcavi are assumed to be occulted by the accretion disk. 
For certain viewing angles there are lines-of-sight which exit 
the grid while in the stellar or disk winds. The winds have 
a much lower density than the cloud material, and therefore 
contribute far less absorption. The model cannot be used 
to calculate the attenuated emission for inclination angles 
smaller than the cavity opening angle due to insufficient ex- 
tent of the computational grid (in reality there is likely to 
be additional absorption from swept-up cloud material for- 
merly in the cavity). 

Once model fluxes have been calculated we use the 
Chandra ACIS-I effective areeE] to convert the synthetic 
spectra into counts space. When calculating the synthetic 
fluxes we assume a distance to the MYSO of 1 kpc and an 
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Figure 3. Snapshots of density (upper panels) and temperature (lower panels) taken from model Rl at (from left to right): t = 1000, 
1200, and 5000 yr. The tick marks on the axis correspond to a distance of 10-*^^ cm. The strip of T ~ 10^ K gas which extends diagonally 
from the wind source is due to frictional heating at the interface between the fast stellar wind and the slower disk wind. 



ISM column (between the observer at Earth and the edge of 
the cloud) of 1.9 x 10^^ cm"^ 

3 RESULTS 

To explore the impact of varying the different parameters 
in the model on both the evolution of the cavity and the 
resulting X-ray emission we have performed an extensive 
parameter space exploration (Table [1]). Each simulation is 
allowed to run until a model time of 2000 yr, except model 
Rl which was allowed to run until a model time of 5000 yr. 
We first describe the cavity evolution and X-ray emission 
features of the fiducial model (Rl), and then examine the 
influence of the various parameters on the results. 

3.1 The standard model 

3AA Dynamical and X-ray properties 

The stan dard model (Rl) mimics the MYSO disk wind 
model of IPrew et al.l (|l998h , with identical wind geometry 
and stellar and disk wind velocities. T he disk wind ma ss- loss 
rate is chosen to match observations (|Felli et al .11 19841 ). The 



stellar wind spee d and mass-loss rate are representative of 
an 08-9V stai0 (|Howarth Prinial [19891 ). The mass infall 
rate for accretion onto the disk, Minfaii, is chosen to be inter- 
mediate between the value ofl.lxlO^'^Moyr"^ derived by 
lAlvarez et aP (|2004al ) for the MYSO Mon R2 IRS 3A, and 
theoretically pre dicted mass infall rates of ^ 10~^ Mp:) yr~^ 
formassive stars (|Baneriee Pudrit j|2007l : iKrumholz et al.l 
l2009l ). The centrifugal radius of the cloud, Tc, is taken to be 
of the order of the observed disk radius for MYSOs (e.g. 
lAlvarez et al.l [2"004al ). As a simple first approximation the 
cavity radius, rcav, is taken to be equal to the centrifugal 
radius of the cloud, and the cavity is assumed to be conical 
(/3=1). 

Fig. [3] shows the spatial distribution of material in the 
simulation; the disk wind lines the cavity wall and separates 

^ The stellar mass adopted in the simulations of 10 is lower 
than that expected for an 08-9V star (~ 30 M©). An increase 
in the mass of the central star will increase the gravitational at- 
traction felt by the infalling cloud material, and thus the infall 
velocity will increase and the density will be decreased. However, 
tests show that this difference does not significantly affect the 
results. 
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Table 1. List of the parameters used in the simulations. /3 is the parameter which defines the curvature of the cavity wall (Eq. [8)), rc 
is the centrifugal radius, uj is the cavity half opening angle, Minfall is the mass infall rate of the molecular cloud, is the transition 
angle between the stellar and disk wind, and is the limiting angle above which there is no wind (i.e. only an accretion flow). For the 
wind terminal velocity, v, mass-loss rate, M, and power, ^, the subscripts "S" and "D" refer to the stellar and disk wind components, 
respectively. = 10 Mq and rcav = 3.0 x 10^^ cm in all models, except models AFGL 2591 and Mon R2 where rcav = 4.5 x 10^^ and 
1.5 X 10-*^^ cm, respectively. 
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the stellar wind from the infalling molecular cloud mate- 
rial. Close to the star there is a small strip of frictionally 
heated wind gas which reaches temperatures of ^ 10^ K and 
therefore emits soft X-rays {E ^ 0.1 keV), although with 
a negligible observed flux when compared to the shocked 
winds. The stellar and disk winds in the simulations are 
supersonic, so that flow incident against the cavity wall re- 
sults in shocked gas bounded by a reverse shock (this is 
best seen in the temperature images in Fig. [3]). The post- 
shock temperature is dictated by the preshock wind speed 
normal to the shock. For vs = 2000 kms~^ , temperatures 
up to T ^ 10® K are obtained, whereas the slower disk 
wind (vd = 400 kms~^) is shocked only up to tempera- 
tures ^ 10^"^ K. At the base of the cavity wall the small an- 
gle of incidence to the normal for the winds, and thus small 
angle of reflection, causes the shocked winds to intersect 
the rotational symmetry axis and close off the unshocked 
stellar wind. A pressure balance between the shocked and 
unshocked gas causes the reverse shock to continually en- 
close the latter. The reverse shock is almost normal to the 
preshock flow near the disk plane and on the symmetry axis 
where the highest post shock temperatures are obtained. Be- 
cause the ram pressure of the inflow/outflow is angle de- 



pendent and the base of the cavity is subject to instabili- 
ties the shape of the reverse shock is often non-spherical. In 
Fig. [3] the reverse shock can be seen in the bottom left of 
the grid, where the enclosed preshock wind has temperature 
T = 10^ K. The position of the reverse shock oscillates and 
its shape changes with time, due to small fluctuations in the 
shape and size of the base of the cavity wall as inflowing ma- 
terial is ablated and incorporated into the outflow, and as 
new inflowing material replenishes it. The size of the reverse 
shock is tied to the postshock pressure. The oscillations have 
a large effect on the observed emission as the postshock den- 
sity, and therefore the emission measure, are dependent on 
the shock position. The postshock stellar wind has velocity 
vectors which are preferentially aligned in the polar direc- 
tion. Small kinks to the shape of the reverse shock cause 
shear velocities in the post-shock flow, and KH instabilities 
are produced which can clearly be seen at t = 1000 yrs. The 
shear layer between the stellar and disk winds provides a 
site for the growth of ^ 10^^ cm amplitude KH instabilities 
on timescales of ~ a few years. By t = 1200 yrs an insta- 
bility of this proportion can be seen driving a clump of disk 
wind material into the path of the stellar wind, which leads 
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to mass- loadi ng of the latt er (for a review of mass- loading 
processes see IPittardllioorl ). 

On larger scales the half opening angle of the cavity, 
cj, decreases until a temporary pressure balance is attained 
between the gas on either side of the cavity wall at t :^ 
1000 yr. A dense layer of compressed cloud material which 
grows with time results from the pile-up of inflowing cloud 
material. At the top of the cavity a small narrowing has 
occured due to the pressure difference across the cavity wall; 
if no winds are blown at all the cavity closes up over the 
star after ^ 500 yrs. Tests performed using the streamline 
equation (Eq. [T]) to calculate the initial cavity shape and 
position did not prevent this occurence. By a simulation 
time of t = 2000 yr, the cavity has widened slightly at the 
base. The inflowing cloud material can still reach the disk 
plane and it is conceivable that accretion onto a disk could 
still be on-going while the winds are interacting with the 
cavity. 

Model Rl was allowed to run for an extended simulation 
time of 5000 yr. At t > 2000 yr the shape of the cavity wall 
continues to change. The base of the cavity wall becomes 
more curved as time goes on, and the cavity wall at the 
top of the grid narrows further (Fig. [3]) . The spectrum at 
t — 5000 yrs has a very similar shape to that at 1000 yrs, 
though with a factor of 2 lower normalization (Fig. [7|). 
Similarly, the Chandra count rate is lower (model Rliong in 
Table [2]) • These differences are unsurprising as the reverse 
shock is now further from the star. 

The X-ray emissivity of the shocked gas varies enor- 
mously with position. The intrinsic X-ray emission comes 
mainly from disturbed cloud material, then the disk wind, 
and finally the stellar wind (Table [2]). The shock driven into 
the cloud is too slow to heat gas up to X-ray emitting tem- 
peratures, but large quantities of cloud material are ablated 
by the outflowing disk wind and mixed into this hotter flow 
(note that the cloud material does not directly mass-load 
the stellar wind - Fig. [3] shows that it is the disk wind which 
does this). This process heats the cloud material to tem- 
peratures where (soft) X-rays are emitted. This process of 
"mass- loading" the outflowing disk wind creates densities 
which are typically two orders of magnitude higher than the 
shocked stellar wind at comparable heights above the disk 
plane. However, the difference between the attenuated lu- 
minosities of the cloud and stellar wind components is very 
small (Lattc = 8 X 10^^ and Lattg = 7 x 10^^ ergs"-^ , respec- 
tively). The explanation is that the stellar wind emission is 
harder and extends to higher energies, and is less affected by 
attenuation. In contrast, the cloud material, which is heated 
to lower temperatures, emits prolifically at low energies and 
has a much higher intrinsic luminosity, but its emission is 
subject to considerable attenuation at £^ < 1 keV. Exam- 
ining the attenuated spectrum at i = 60° (Fig. 2]), we find 
that the stellar wind material dominates the emission at 
^ > 4 keV, with comparable contributions from the winds 
and cloud material below this energy. 

The observed model X-ray emission for i = 60° shows 
variability on timescales of order a few hundred years 
(Fig. [5|). The total emission initially rises after the stellar 
and disk winds are switched on and a reverse shock is al- 
lowed to develop. There are then quasi-periodic oscillations, 
with noticeable dips at t 500, 1000, 1300, and 1700 yrs. 
The fluctuation in the total emission approximately follows 




Energy (keV) 

Figure 4. Attenuated 0.5-10 keV X-ray spectra for model Rl 
at various inclination angles. From top to bottom: i = 30,45, 60, 
and 90°. The total (red), stellar wind (green), disk wind (blue), 
and heated cloud material (pink) emission are shown. All spectra 
were calculated at t = 1000 yr and have been convolved with the 
Chandra effective area. Note the difference in scale in the top 
panel. 
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Figure 5. Attenuated 0.1-10 keV lightcurve for model Rl with 
i = 60°. The hghtcurves of the separate emission components 
(stellar wind, disk wind, cloud) are also plotted. 



the behaviour of the emission from the stellar wind. Some 
features are due to the disk wind luminosity, which is some- 
times the dominant emitter. Fluctuations in both of these 
emission components is strongly related to the position of 
the reverse shock. At t — 500 yr the distance of the re- 
verse shock from the star increases, and a reduction of the 
stellar wind luminosity occurs. Similarly, the maximum to- 
tal luminosity occurs at t = 1200 yr, in conjunction with 
a minimum in the distance of the reverse shock from the 
star and a shape which results in more stellar wind ma- 
terial being shocked to the highest temperatures (Fig. [3]). 
The variability in the observed emission is a general feature 
for all the simulations and complicates comparisons between 
the different model results. In simulation Rl the cavity has 
reached an approximate steady state and therefore the emis- 
sion averaged over 500 or 1000 yrs should remain roughly 
constant. At t = 1200 yr a clump of disk wind is driven into 
the stellar wind by turbulent mixing (Fig. [3|). Such mix- 
ing boosts the observed emission of both winds by heating 
the cooler disk wind and increasing the emission measure 
of the stellar wind. However, the bulk of the emission still 
originates from close to the reverse shock. The temperature 
of the stellar wind gas is also increased as it impacts the 
clumps and additional shocks are produced. At t > 2000 yr 
the attenuated luminosity appears to steadily decrease and 
the amplitude of the fluctuations seem reduced. This sug- 
gests that the variability at t < 2000 yr may be related to 
an initial evolutionary stage and thus to the initial condi- 
tions. Whether the variability to the observed emission will 
continue on longer timescales (i.e 10^ yrs) is not clear as 
it depends on the postshock gas pressure. The luminosity 
of the cloud material follows the morphology of the winds, 
which suggests the magnitude and variability of this emis- 
sion is related to the activities of the winds. 

The visual extinction to a deeply embedded star relates 
to the column density along that line-of-sight. To allow com- 
parison between our calculated column densities and obser- 
vationally determined visual exctinctions, Av's, we note the 
simp le convers ion factor — A^h/(1-9 x 10^^ cm~^) mag 
(e.g. ICoxll2000l V One method of observationally determining 
the optical depth to the MYS O is to use the 9.7 /im Sili- 
cate feature ([Crapsi et ahllioosl v For the simulations, a good 
approximation is the column density to the star, A^H-star 



(i.e. the grid origin). Comparing this value to the emission 
weig hted column (EWC), iVEwc = SA^Hci^intTc/SLintTc 
(e.g. IParkin Pittardll2008l ) , where the index c indicates the 
column and luminosity from each cell on the hydrodynamic 
grid and the summation is over all cells, shows that the col- 
umn density to regions with large intrinsic X-ray emission is 
significantly higher than that to the star (Table [2]). This is 
mainly due to lines-of-sight to the X-ray emitting plasma in 
the receding lobe having longer path lengths and/or passing 
through denser material. For an inclination angle oil = 60°, 
the majority of the absorbing column is through the first 
2000 au of cloud gas beyond the cavity wall with a large 
fraction due to the layer of compressed cloud material which 
lines the cavity wall. 

3.1.2 Inclination angle dependence of the observed X-ray 
emission 

The attenuation of the emission is viewing angle depen- 
dent. When the inclination angle \s i = 30° the majority 
of lines-of-sight to X-ray emitting plasma avoid the dens- 
est cloud material close to the equatorial plane. This re- 
sults in a lower A^H-star and significantly brighter spectrum 
at < 4 keV than at i = 60° (Fig. At E > A keV 
the harder stellar wind component dominates the spectrum. 
The disk wind and heated cloud material contribute the 
majority of the X-ray emission, with little difference be- 
tween the two (Latto — 1-24 x lO^^ergs"^ and Lattc — 
1.61 X lO^^ergs"^ , respectively). The total attenuated lu- 
minosity (LattT — 2-88 X 10^^ ergs~^ ), and the integrated 
0.1-10 keV count rate (77 :^ 30 ks~^) are two orders of mag- 
nitude higher than observationally detected count rates. 

When i — 45°, the majority of the X-ray emitting 
plasma is viewed through the cloud material. The column 
densities have increased relative to i = 30° as now the 
line-of-sight to the star and to the regions with highest in- 
trinsic emission must pass through sections of higher den- 
sity, compressed cloud material. Examining the spectrum 
(Fig. |4]) the turnover energy is now at ^ :^ 1.5 keV. This 
is an increase relative to the i — 30° case which illus- 
trates the higher attenuation to the X-ray emission from 
disk wind and cloud material. At energies E ^ A keV the 
spectrum has not changed significantly. Lattx ^ind 77 are now 
more in line with X-ray observation s of massive star form- 
ing regions (e.g. iKohno et al. I I2OO2I : IPreibisch et al. 1 12OO2I : 
iGiardino etaD l2004l ) . Increasing the inclination angle fur- 
ther to i = 60° results in a slight rise in A^n-star, A^ewc, 
and the turnover energy of the spectrum. Similarly, LattT 
and Tf decrease. 

The maximum column densities are attained when 
viewing the system edge-on (i.e. i — 90°). It is interesting 
to note that the column to the star is now greater than the 
EWC, which indicates that for this viewing angle the column 
to the spatially extended (few thousand au) distribution of 
X-ray emitting plasma is lower than that to the star. Un- 
like the lower inclination angle calculations, emission below 
E ^ 2 keV is now observed from both lobes of the cavity. 
Also, the stellar wind material, which has the lowest intrinsic 
luminosity, now provides the largest contribution to LattT- 
The spectrum above E ^ 2 keV remains roughly similar. 
We note that for inclination angles where the line-of-sight 
to the X-ray emitting plasma passes through the cloud ma- 
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Figure 6. Broadband images of the attenuated X-ray emission from model Rl at t = 1000 yr, an inclination angle of i = 6i 
of 1 kpc, and in the energy bands 1-2 (left panel), 2-4 (middle panel) and 4-10 keV (right panel). 



terial (i.e. i ^ 45° when uj 30°), emission above E ^ 4 keV 
is mainly contributed by the shocked stellar wind. 

The inclination angle dependence of the observed emis- 
sion shows that there is a narrow range of angles where the 
count rate exceeds the Chandra X-ray detections. 

3.1.3 Synthetic X-ray images 

The raytracing code can produce synthetic broadband X- 
ray images which allow us to identify where different energy 
X-rays preferentially originate from. For simulation Rl at 
i = 60°, we see that the observable X-ray emission in the 
1-2, 2-4, and 4-10 keV bands originates from similar regions 
of the cavity (Fig. [H]). A peak in the intensity in the three 
bands originates near the reverse shock, and is mainly gen- 
erated by shocked stellar and disk wind material (Fig. |4]). 
Emission extending in the polar direction traces the bound- 
ary between the hot shocked stellar wind, and the cooler, 
higher density disk wind. 

Examining Fig. [3] shows that the region of the cavity 
close to the pole is filled by hot shocked stellar wind. This 
gas emits the majority of the emission above £^ ^ 4 keV 
(Fig. 2]). The tower of emission visible in the 4-10 keV image 
in Fig. [6] therefore traces the stellar wind, with the shocked 
cloud emission in this energy band originating from the base 
of the cavity wall. The peak in intensity in all three bands 
occurs at the waist of the reverse shock, where there is a 
combined maximum in temperature and density. Examining 
the spectrum in the E = 2 — 4 keV energy range shows that 
this emission is from shocked disk wind and cloud material, 
which suggests it originates from the T ^ 10^"^ K gas at 
z = r 0.5 X 10^^ cm. Interestingly, both the preceding 
and receding lobes are observable in the 2-4 and 4-10 keV 
bands, which signifies the lessening influence of attenuation 
as the photon energy increases. 

Importantly the spatial extents of the emission in all 
three energy bands (1-2, 2-4, and 4-10 keV) are just below 
the resolution limit of Chandra 0.5^^) and so this model is 
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consistent with the unresolved nature of real MYSOs. How- 
ever, this finding is not true for all of our models which 
allows us to place constraints on some of the key model pa- 
rameters. 

3.2 Variation with mass inflow/outflow rates 

The velocities of the inflowing and outflowing gas in the 
simulations are supersonic. As such, the morphology of the 
evolving cavity is strongly dependent on the ratio of the in- 
flow and outflow ram pressure, which is directly proportional 
to the mass flow rates. There are a multitude of possible 
models with identical ram pressure ratios, but producing 
vastly different observable X-ray emission characteristics. 
Since optically thin bremsstrahlung emission has a density 
squared dependence (the emissivity oc n^), an increase in the 
mass-loss rates can have a considerable impact on the ob- 
served flux. However, an increase in the emission from cloud 
material due to an increase in Minfaii may be countered by 
an accompanying increase in attenuation. The degree of ab- 
sorption is dependent on the column density, which is largely 
due to the cold cloud material, and therefore is directly pro- 
portional to the mass infall rate (A^h oc Minfaii)- 

Decreasing the mass infall rate reduces the ram pres- 
sure of the infall. In model R2 the mass infall rate is halved 
to Minfaii = 1 X 10~^ Mq yr~^ . This results in less resistance 
to the outflows and subsequent expansion of the cavity. The 
cavity wall is hence wider at the base at an equivalent time. 
Also, the reverse shock is positioned at a greater distance 
from the star, and the enclosed region appears more elon- 
gated in the 2;-direction. Apart from the afore- mentioned 
differences we note that simulations R2 and Rl show sim- 
ilar evolution. As expected, A^H-star and A^ewc decrease 
by roughly a factor of two, consistent with the decrease in 
Minfaii- There is a slight increase in LintT but LattT and 77 
increase by factors of 2 and 3.5 respectively. Compared to 
model Rl, the spectrum for model R2 (Fig.[7|) shows higher 
flux at £^ < 4 keV and a softer spectrum above this en- 
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Figure 7. Total attenuated 0.5-10 keV spectra for simulations. 
The top panel shows Rl, Rliong^ R2, R3, R5, the upper middle 
panel shows Rl, R7, R8, R9, RIO, the lower middle panel shows 
Rl, R12, R13, and the lower panel shows Rl, R15, R17, R18. All 
spectra were calculated at a time of t = 1000 yr (except models 
Rliong which was at t = 5000 yr, and R3 which was at t =500 
yrs), using an inclination angle, i = 60°, and have been convolved 
with the Chandra effective area. Note the difference in scale in 
the upper middle panel. 




Figure 8. Density snapshot taken from model R4 at t = 2000 yr. 
The tick marks on the axis correspond to a distance of 10 -"^^ cm. 
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Figure 9. Density snapshots taken from model R7 at t =1700 
(left) and 2000 yrs (right). The tick marks on the axis correspond 
to a distance of 10-*^^ cm. 



ergy which is linked to the differences in the position of the 
reverse shock. 

If the mass infall rate is increased the ram pressure of 
the infalling material increases in direct proportion and the 
outflowing winds encounter more resistance to their expan- 
sion. In model R3 the reverse shock is pushed closer to the 
star by the increased mass infall rate, and the cavity now 
closes up over the star. At t = 500 yr the infall has begun to 
obstruct some of the outflow, and by t = 900 yrs the outflow 
has been fully constricted. Table [2] contains X-ray luminosi- 
ties at t = 500 yr, before the cavity closes up (once closed 
the X-ray emission in our model becomes zero). Examining 
Table [21 we find that A^H-star and A^ewc have increased, 
and the Latt's and ?7's have decreased for R3 relative to Rl. 
Note that Lintc, increases considerably with Minfaii, as the 
density, and thus the emission measure of this material, in- 
creases. However, the increasing attenuation with Minfaii re- 
moves more and more of the soft X-ray emission, so that the 
observed emission becomes increasingly hard (see Fig. [7]). 
Model R3 shows that there is a delicate hydrodynamic bal- 
ance between the infall and outflow, and large contrasts be- 
tween these parameters can result in vastly different cavity 
evolution. 

Confinement of a wind was dis cussed in the context 
of T-Tauri stars by IChevalierl (Il983h . where the infall ram 
pressure can suppress the expansion of a wind into the sur- 
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rounding envelope. We have a similar situation, albeit more 
complicated by the anisotropic outflow. The wind at the 
base of the cavity may have a higher ram pressure than 
the stellar wind towards the pole of the cavity (as in model 
R3), which may ultimately lead to confinement at the poles 
but expansion at the base. Note also that the simulations 
we have performed only consider hydrodynamics, whereas a 
number of physical processes may in reality prevent a cavity 
from closing up. Magnetic fields may play an important role 
in controlling the dynam ics of the infalling envelope (e.g. 
IBaneriee Pudrita 12007), or the radiatio n pressure may 
halt the inflow (i.e. iKrumholz et al.ll2005h . However, it is 
clear that within the context of our model, the parameters 
for a slowly expanding cavity lie within a narrow range. 

In model R4 the mass-loss rate of the stellar wind is 
increased to 1 x 10~^ M© yr~^ . The boundary between the 
inflowing and outflowing gas is likely to be subject to dy- 
namical shear instabilities. Unlike in Rl, such instabilities 
can be seen forming in the lower part of the cavity wall after 
t ^ 500 yr, with prominent structure visible at t = 2000 yrs 
(Fig. [8]). A common feature in the simulations where these 
instabilities are noticeable (R4, R5, R9, R14, R15, and R18) 
is that the disk wind is compressed into a thin layer along 
the cavity wall. The timescale for KH instabilities to grow 
is oc {pi + p2)l ^ pip2 , and the enhanced disk wind density 
along the cavity wall results in an increased growth rate for 
these instabilities. In addition to this, in the afore-mentioned 
simulations the cavity wall becomes curved at the base due 
to the action of the winds (except in R18 where the cavity 
is initially curved). As such, the shocked winds have a small 
grazing angle against the base of the cavity wall and the ve- 
locity difference across the shear boundary is higher. These 
instabilities cause the initially smooth cavity wall to become 
deformed, and in the more severe cases these deformations 
grow into jagged features. As these features protrude into 
the winds clumps of cloud material are stripped away and 
then ablated. In model Rl, we see disk wind material mixing 
into stellar material (see Fig. [3]), while in model R4 (Fig. [8]) 
we see cloud material mixing into disk wind material. Such 
mixing can cause the disk wind to become more confined to 
the lower regions of the simulation domain, and no longer 
border the entire cavity wall. Also of note in model R4 is 
that the wind-driven expansion of the cavity opposes the in- 
fall and it is conceivable that at t > 2000 yr the flow of cloud 
material onto the disk plane will be halted and the envelope 
will be completely destroyed by the winds. The increase in 
Ms results in an increase in the Latt's and r\. 

In model R5 the ratio of the inflow/outflow ram pres- 
sures is the same as in R4, but Minfaii, Ms, and Md are 
all ten times lower. Nevertheless, the evolution of the cavity 
and the position and shape of the reverse shock are almost 
identical between R4 and R5, as one would expect. The in- 
trinsic emission decreases by a factor of 100, as expected 
[EM oc M^). However, r\ only decreases by a factor of 2, 
since the attenuation is also vastly reduced. Comparing the 
spectra (Fig.[7| we see that in the 4-10 keV band the spectral 
slope is almost identical, yet for R5 the count rate in this 
band is roughly a factor of 100 lower. Therefore, although 
count rates consistent with observations can still be attained 
when the mass flow rates are scaled down, the much softer 
observed spectrum places constraints on the outflow, and 
consequently the inflow, parameters. 




Figure 10. Density snapshots taken from model R8 at f = 
600 (left) and 900 yrs (right). The tick marks on the axis cor- 
respond to a distance of IQ-*^^ cm. 



In R6 the mass-loss rate of the disk wind is increased by 
a factor of 3 from Rl to Md = 3 x 10~^ M© yr"^ . If Mnfaii 
is kept at 2 X 10~^ M© yr~^ the disk wind drives a hole 
into the base of the cavity wall. To prevent this occurence 
Minfaii is increased to 1 x 10~^ M© yr~^ . The resulting cavity 
evolution is somewhat similar to that of Rl. However, the 
fact that Minfaii has been increased to stabilise the cavity 
evolution brings an accompanying increase in the column 
densities so that the attenuated luminosity is actually lower 
than that from model Rl. Compared to model R3 (which 
has identical parameters except for a lower value of Md) 
we see that the higher value of Md prevents the infall from 
totally overwhelming the outflow. 

Reducing the mass-loss rate of the stellar wind (R7) 
leads to major differences in the cavity evolution. The pres- 
sure applied to the cavity wall by the winds as a function 
of polar angle now differs more than in model Rl (there 
is now a factor of 20 variation in Mq). Whereas the stellar 
wind previously had sufficient ram pressure to hold back the 
infalling cloud material, the cavity now begins to close up 
above the star. At t ~ 1500 yrs the disk wind begins to 
cut into the base of the cavity wall, and the bubble formed 
continues to increase in size for the duration of the simula- 
tion (Fig. [9]) . By t = 2000 yrs a limb of cloud material is 
extending towards the star. Further analysis shows that by 
t = 2300 yrs this material breaks away from the cavity wall 
and obstructs the outflowing winds close to the star. 

3.3 Variation with wind speeds 

The preshock velocity dictates the postshock gas tempera- 
ture (T oc v'^). Hotter gas emits higher energy X-rays, and 
an increase or reduction in the wind speeds will affect the 
hardness of the resulting spectrum. 

In model R8 the stellar wind velocity is increased to 
— 3000 kms~^. The shocked stellar wind now reaches 
T > 10^ K. The equatorial flow appears to be more powerful 
and drives a small, short-lived hole into the base of the cavity 
wall (Fig. I10|) . This is unexpected as it is the ram pressure 
of the stellar wind which has been increased. However, it 
seems that the total pressure near the equatorial plane is 
enhanced by an increase in the thermal pressure adjacent 
to the reverse shock. Subsequently, the cloud material then 
collapses (see Fig. llOp and fills the hole, after which the 
cavity settles into a steady state similar to that of model 
Rl. The higher ram pressure of the stellar wind increases the 
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Table 2. Column densities (columns 3 and 4), intrinsic (int) and attenuated (att) 0.1-10 keV X-ray luminosities (columns 5-8 and 9-12, 
respectively), and integrated count rates (column 13) from raytraced emission calculations. The stellar wind (S), disc wind (D), and 
cloud material (C) contributions to the total (T) luminosity are shown. The integrated 0.1-10 keV count rate, 77, was calculated using 
the Chandra effective area (see ^ I2.5|) . Luminosities are in 10^-*^ ergs"-*^ , and count rates are in ks"-*^. A^H-star and A^ewc a-^e the column 
density measured to the star and the emission weighted column respectively, both of which are in 10^^ cm~^. The raytracing calculations 
were performed with an inclination angle of i = 60°, except RI30, RI45, and RI90 which were calculated using inclination angles of 0, 
30, 45, and 90° respectively. All calculations were performed on simulations at t = 1000 yr, except Rliong and R3 which were calculated 
at 5000 and 500 yrs, respectively. 
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Lints 


LintD 


Lintc 


LattT 


Latts 


Latt^ 


Lattc 


J] 


Rl 




2.47 


6.09 


370 


0.500 


43.1 


330 


0.019 


0.007 


0.004 


0.008 


0.11 


RI30 




0.44 


5.42 


370 


0.500 


43.1 


330 


2.88 


0.035 


1.24 


1.61 


30 


RI45 




3.11 


6.21 


370 


0.500 


43.1 


330 


0.040 


0.008 


0.012 


0.020 


0.35 


RI90 




9.66 


7.03 


370 


0.500 


43.1 


330 


0.012 


0.006 


0.002 


0.004 


0.05 


Rllong 


longer run 


2.40 


5.54 


186 


0.221 


26.1 


157 


0.007 


0.002 


0.004 


0.001 


0.05 


R2 


illidii 


0.99 


3.19 


436 


0.313 


32.2 


404 


0.040 


0.005 


0.021 


0.014 


0.37 


R3 


Minfall T 


16.7 


46.4 


3070 


0.596 


27.7 


3050 


0.014 


0.002 


0.005 


0.007 


0.02 


R4 


Ms T 


1.92 


3.46 


220 


1.39 


33.4 


185 


0.071 


0.043 


0.015 


0.012 


0.41 


R5 


R4 with M's i 


0.25 


0.37 


2.48 


0.11 


0.22 


2.16 


0.016 


0.002 


0.006 


0.008 


0.23 


R6 


Minfall T, Md T 


16.9 


46.6 


6120 


5.27 


395 


5720 


0.008 


0.007 


7 X 10-4 


5 X 10-4 


0.02 


R7 


Ms i 


2.47 


6.39 


390 


0.470 


44.1 


345 


0.009 


0.005 


0.002 


0.002 


0.05 


R8 


vs T 


2.48 


7.38 


569 


0.142 


33.8 


535 


0.049 


0.005 


0.003 


0.011 


0.26 


R9 


i, Ms T 


2.10 


4.40 


215 


4.49 


36.8 


174 


0.042 


0.027 


0.008 


0.007 


0.28 


RIO 


B star 


0.38 


0.96 


13.5 


0.066 


2.26 


11.2 


0.004 


4 X 10-4 


7 X 10-4 


0.003 


0.06 


Rll 


Os = 70° 


2.83 


7.42 


438 


1.64 


31.5 


405 


0.019 


0.013 


0.003 


0.003 


0.08 


R12 


Os = 85°, On = 89° 


3.00 


9.15 


526 


2.41 


6.93 


516 


0.030 


0.017 


0.002 


0.012 


0.12 


R13 


no disk wind 


2.93 


9.52 


456 


2.22 




454 


0.060 


0.032 




0.029 


0.20 


R14 


cc; = 5° 


2.41 


6.66 


318 


0.253 


27.8 


290 


0.051 


0.011 


0.016 


0.024 


0.30 


R15 


LJ = 10° 


2.82 


7.22 


483 


0.577 


54.5 


4.28 


0.009 


0.006 


0.002 


0.002 


0.05 


R16 


u; = 45° 


2.14 


6.40 


219 


0.219 


23.8 


195 


0.018 


0.004 


0.012 


0.002 


0.12 


R17 


Lj = 60° 


1.21 


3.36 


27.5 


0.125 


11.9 


15.4 


0.003 


2 X 10-5 


3 X 10-4 


0.003 


0.05 


R18 


f3 = 2 


1.36 


3.40 


341 


1.11 


30.3 


310 


0.099 


0.036 


0.015 


0.048 


0.58 


R19 




2.86 


7.42 


754 


0.378 


45.2 


708 


0.025 


0.006 


0.009 


0.010 


0.14 


R20 


ret 


2.39 


6.21 


148 


0.158 


11.7 


136 


0.079 


0.008 


0.059 


0.012 


0.39 



maximum distance of the reverse shock from the star, and its 
amplitude of oscillation. Increasing the stellar wind velocity 
causes the spectral hardness (Fig. O and the Chandra count 
rate, 77, to increase. Interestingly, the X-ray emission from 
the disk wind and cloud material is now dominant at all 
energies. Whereas Latts Lattc were approximately equal 
in model Rl, Lattc twice as high as Latts model R8. 
This is likely a result of increased turbulent heating by the 
stellar wind. A decrease in Lints occurs because the reverse 
shock is further from the star, leading to a lower post shock 
density, and thus a lower emission measure. 

Recent theoretic al models of stellar evolution by 

iHosokawa Omukail (l2009l ) predict that for high accretion 
rates (~ lO~^M0yr~^) the protostellar radius becomes 
very large. This results in a lower surface gravity for the 
protostar and consequently the stellar wind will be more 
like that of a supergiant (i.e. vs ^ 1000 kms~^) than 
a main sequence 0-star. Model R9 explores this scenario, 
where Ms = Md = lO~^M0yr"^ and ^;s 1000 kms"^ . 
With a greater ram pressure in the stellar wind compared 
to model Rl, the cavity does not initially narrow at the top 



of the grid. The postshock stellar wind has a lower temper- 
ature (T ^ 10^"^ K) because of the lower preshock velocity. 
The position of the reverse shock is much further from the 
star [z ^ 4: X 10^^ cm), and this distance increases as the 
simulation progresses. The column densities are lower (Ta- 
ble [2]) , likely due to the greater height of the reverse shock 
(and thus the postshock gas) above the disk plane, and thus 
the reduced path length through cloud material for lines-of- 
sight. The spectrum is very different from that of Rl and 
is clearly a lot softer (Fig. [T]). The stellar wind component 
now dominates the observed emission at L^ 1 keV as the 
denser stellar wind increases the emission measure. The syn- 
thetic X-ray images show the approaching lobe to be fully 
illuminated by X-ray emission, unlike the narrow tower seen 
for Rl in Fig.[6l with the intensity peaks in the 1-2, 2-4 and 
4-10 keV bands being > I" in extent (see Fig. [TT]) . 

Comparing models R4 and R9 gives insight into the ef- 
fects of reducing the stellar wind speed, vs, with all other 
parameters kept the same. In a similar respect to models Rl 
and R8, the reduction in vs causes the reverse shock to be 
closer to the star, and the postshock stellar wind tempera- 
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Figure 11. Broadband image of the attenuated X-ray emission 
for model R9 at an inclination angle of i = 60°, a distance of 1 
kpc, and in the 4-10 keV energy band. 

ture to be lower (~ 10^"^ K). Column densities are slightly 
higher for model R9. The Lint's are very similar, yet the 
Latt's and 77 are slightly lower. Comparison of the spectra 
shows that for model R9 the flux level is slightly lower and 
the spectral slope at L^ > 4 keV is slightly softer than model 
R4. 

A number of known objects possessing outflow cavities 
where X-ray emission has been detected have been inferred 
as early B-type stars. These stars will have lower termi- 
nal wind speeds and mass-loss rates than the 0-type star 
used as our fiducial model. In model RIO we use stellar 
wind parameters corresponding to a BIV star (Prinja 1989). 
Consequently, the mass infall rate must be scaled down to 
produce a similar ram pressure balance to that of model 
Rl. The evolution of the cavity is qualitatively similar to 
model Rl. However, the different ram pressures of the in- 
flow and outflow produces more vigourous ablation of the 
base of the cavity wall than in model Rl. The reduced stel- 
lar wind velocity results in a lower maximum postshock gas 
temperature, which in-turn reduces the amount of hard X- 
ray emission. The column densities are significantly reduced 
compared to model Rl, and are in agreement with observa- 
tion ally determined value s for embedded early B-type stars 
(e.g. IPreibisch eralll2002l ) . The spectrum is much softer now 
(Fig. El). 

3.4 Variation with wind geometry 

Varying the wind geometry affects the contributions to the 
total emission from the different components, and changes 
the cavity evolution. For instance, increasing the opening 
angle of the stellar wind provides more stellar wind material, 
and less disk wind material. If the ram pressure of the stellar 
wind is much greater than that from the disk wind the latter 
may be channelled into an equatorial outflow, as found in 
model R8. To explore the consequences of different wind 
geometries we have calculated models Rll, R12, and R13 
(see Table p. 



As the opening angle of the stellar wind is increased 
(c.f. models Rl, Rll, and R12), and that of the disk wind 
decreased, the thickness of the disk wind layer along the 
cavity wall becomes thinner due to the reduced amount of 
disk wind being injected. The evolution of the cavity remains 
similar to that of Rl, although the cavity wall is now more 
curved towards the base. The base of the cavity wall is closer 
to the star because the ram pressure of the outflow close to 
the disk plane is now reduced. The aspect ratio of the region 
of unshocked winds increases as the reverse shock moves to 
smaller r. The value of 77 does not vary much from that 
calculated for Rl though the observed spectrum is harder 
(see the spectrum from model R12 in Fig. [7|). 

To explore the effect of having a spherical outflow we 
have constructed model R13, which is equivalent to Rl ex- 
cept with no disk wind. We find that compared to model 
Rl there is a definite hardening with more X-ray emission 
at £^ > 4 keV (Fig. O, which is expected as there is now 
more high speed material in the outflow. There is also an 
enhancement in Lattc due to heating of this gas to hotter X- 
ray emitting temperatures by direct contact with, and thus 
greater turbulent heating from, the stellar wind. The value 
of Tf remains in agreement with observationally detected val- 
ues for MYSOs. Interestingly, this shows that X-ray emission 
can be generated from the wind-cavity interaction in the ab- 
sence of a disk wind, which means this mechanism could be 
applicable to more evolved protostars as well. 



3.5 Variation with cavity opening angle 

The angle that the outflow makes with the cavity wall af- 
fects the rate at which kinetic energy in the preshock gas 
is transferred into thermal energy in the shocked gas. Also, 
the pressure of the shocked gas is affected by the degree of 
confinement provided by the cavity wall. Keeping all other 
parameters the same, a narrow cavity leads to higher gas 
pressures than a wider cavity. This has implications for the 
wind-driven expansion of the cavity and the position and 
shape of the reverse shock. To explore these effects further 
models R14-R17 have been performed. 

When ix) is reduced relative to Rl (as in models R14 
and R15) the cavity is initially narrower, and the postshock 
thermal pressure higher. This results in a rapid expansion 
of the cavity in order to balance the pressure across the cav- 
ity wall. In model R14 the cavity seems to have relaxed by 
t = 1000 yrs, whereas in model R15 expansion still appears 
to be ongoing at this time. By the end of both simulations 
the cavity has a roughly curved shape. The rapid expansion 
of the cavity wall away from the star causes turbulent mix- 
ing and fluctuations in the shape and position of the reverse 
shock. Instabilities form in the cavity wall as the less dense 
outflow collides with the more dense inflow, causing a combi- 
nation of KH, Ray leigh- Taylor, and Richtmeyer-Meshkov in- 
stabilities. The column densities relative to those from model 
Rl are higher because of the longer path lengths for lines-of- 
sight through cloud material and because the initial cavity 
is cut into the cloud closer to the star. The conclusion that 
the values of 77 and LattT are higher for smaller initial cavity 
opening angles appears to only be true for model R14. How- 
ever, the lower 77 at t = 1000 yr resides at a minimum in the 
lightcurve for model R15. The spectrum remains similar to 
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Figure 12. Density snapshot taken from model R17 at t = 
1000 yr. The tick marks on the axis correspond to a distance 
of lO^^ cm. 
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Figure 13. Density snapshots taken from model R18 dX t = 1400 
(left) and 2000 yrs (right). The tick marks on the axis correspond 
to a distance of 10^^ cm. 



that of Rl at £^ < 2 keV, and shows more emission above 
this energy (Fig. [7|). 

In models R16 and R17 the cavity half opening angle 
is increased compared to model Rl. For uj = 45° (R16), 
a reverse shock is still generated close to the star which 
encloses the unshocked winds. Interestingly, when co = 60° 
(R17) this no longer happens (see Fig. ll2|) . This is the result 
of a shallow grazing angle for the wind-cavity collision. The 
shocked gas now only resides along the cavity wall, and the 
maximum temperatures reached are T - 10^-^ K. The disk 
wind forms a layer along the cavity wall. On the interior 
side of the cavity it is bordered by the shocked stellar wind. 
Column densities and the intrinsic and attenuated emission 
decrease with increasing uj. When there is no reverse shock 
(close to the star) 77 drops significantly and the spectrum 
becomes much softer (see the spectrum from model R17 in 

Fig. ED. 

3.6 Variation with cavity curve 

The shape of the cavity may not initially be conical, as has 
been assumed in simulations R1-R17 and R19-20. In fact, 
some observations of cavities around MYSOs show them to 
be more parabolic in shape. It is unclear whether an ini- 
tially conical outflow cavity evolves to become more curved. 
Such evolution could be driven by the action of the outflow- 
ing winds (as in Rl), or by changes to the mechanism that 
drives the outflow causing the resulting cavity to be curved 
at the base. The initial curvature of the cavity wall is con- 
trolled by the exponent f3 in Eq.[8] In model R18 we adopted 
/3 = 2, i.e. a parabolic cavity. A reverse shock is still formed, 
although it occurs at a greater distance from the star, and 
with a different shape than that seen in model Rl fFig. I13|) . 
The reverse shock can be divided into two parts. The up- 
per part (which spans polar angles < ^ < 60°) is very 
clearly normal to the stellar wind, which leads to post shock 
gas at very high (^ 10^ K) temperatures and hard and lu- 
minous X-ray emission (Fig.fT]). The lower part, in contrast, 
is very oblique to the disk wind, and leads to much cooler 
postshock temperatures. For the duration of the simulation 
infalling cloud material can still reach the disk plane. How- 
ever, by t = 2000 yr this accretion flow has halted due to 
the pocket of hot rarefied disk wind gas which sweeps up 
the cloud material and drives it into the equatorial outflow 
cavity fFig. ll3|) . Also noticeable by the end of the simulation 
is that the disk wind is confined to the equatorial outflow 



bubble, and the larger cavity towards the pole is filled with 
turbulent shocked stellar wind gas. 

In conclusion, introducing curvature into the cavity wall 
at the start of the simulation does not prevent integrated 
count rates in approximate agreement with observations 
from being produced and the spectrum appears largely sim- 
ilar, though a factor of 3 or so brighter (Fig. [7|). Curvature 
produced by a previous outflow may be erased as the winds 
carve into the cavity wall. 



3.7 Variation with centrifugal radius 

The spherically averaged density goes as p oc r~^^^ within 
Tc, and p oc r~^^^ at distances greater than Tc. Therefore, 
changing affects the density distribution in the cloud. 

In model R19, is decreased to 0.8 x 10^^ cm, and as 
such the density in the cloud material should be lower than 
that of model Rl at equivalent distances from the star. There 
is however little noticeable difference between the cavity evo- 
lution in this simulation and Rl. There is a minor increase 
in A^H-star which is consistent with the cavity wall being 
slightly closer to the star for model R19. The increase in 
Newc is a result of the spatial distribution of the emission 
and the position of the reverse shock. The reverse shock is 
now closer to the star than in model Rl, and as such there 
is more emission from the disk wind close to the base of the 
cavity wall. Lines-of-sight to this disk wind gas pass through 
the dense cloud gas close to the axis, resulting in higher col- 
umn densities. However, in spite of this 77 is higher because 
of the enhancement to the emission measure brought about 
by the position of the reverse shock. 

The value of rc is increased to 6.0 x 10^^ cm in simula- 
tion R20. The winds now excavate the base of the cavity wall 
more. This excavation proceeds in an unsteady manner, and 
results in more variation in the degree of mixing between 
the winds and the position of the reverse shock throughout 
the simulation than seen in model Rl. These differences in 
the dynamics result in increased variability of the resulting 
emission; for model Rl the largest amplitude oscillations in 
the attenuated luminosities was a factor of 2 (Fig. [5]), 
whereas for model R20 there are oscillations of a factor of 
^ 4. The spectrum at t = 1000 yrs represents the median 
for model R20. 
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Table 3. Model parameters (see also Table [T} an d results for the 
candidate objects. The dist ances are taken fro m FSchneider et al.l 
(l2006l ). IStaude et al.l (Hosi) andlRacind (196^ and the inclina- 
tion a ngle s are from van der Tak et alT (l2006l i. ISolf Carsentvl 
(Il982l ) and lAlvarez et al.l (l2004al ) for AFGL 2591, S106 IRS4, and 
Mon R2 IRS3 A, respectively. A^H-star and Aewc are in units of 
1023 cm-2. 



Model 


D 


i 


A^H-star 


Aewc 


V 




(kpc) 


n 






(ks-i) 


S106 


0.6 


75 


0.4 


0.8 


0.23 


Mon R2 


0.83 


45 


0.6 


1.1 


0.12 


AFGL 2591 


1.7 


32 


2.0 


6.9 


0.09 



4 CANDIDATE OBJECTS 

We compare our model against X-ray observations of three 
objects: S106 IRS4, Mon R2 IRS3 A, and AFGL 2591. Be- 
cause of the low count rates spectral fits to these objects 
are relatively meaningless. Therefore we restrict our com- 
parison to the integrated count rates and the generation of 
X-rays above 2 keV. We also examine whether the model 
flux falls below the detection limit (and is spatially unre- 
solvable) when reasonable model parameters are used. We 
place the constraint on the model parameters that the cavity 
must be reasonably steady (i.e. the infall does not constrict 
the outflow at any point during the simulation) and the vi- 
sual extinction agrees with observation. Therefore, the main 
free parameters are Ms and Minfaii (as disk wind parameters 
were based on previous observational estimates). 

4.1 S106 IRS4 



S106, at a distance of 600±100 pc (IStaude et al.lll982l ). is a 
massive star-forming region known for an extended bipolar 
HII region, whic h is illuminated b y the 15 massive 
star S106 IRS4 (iFelh et al.|[T983 ). ISchneider et all (l2002l ) 
find the observed morphology and kinematics in ^^CO and 
^^CO 2^1 to be consistent with a cavity created by the 
radiation and ionized wind from S106 IRS4 sweeping up 
material from the cavity wall, marked by the two lobes of 
the HII region and at the extreme ends of the flow. They 
also find that the double-peak structure of the cloud breaks 
down and the emission merges into a more diffuse extended 
plateau at S106 IRS4, where the molecular emission traces 
the red-shifted component of the stellar wind hitting the 
backside of th e cavity walls. Near-IR speckle observations by 
I Alvarez et al.l (|2004bl ) show o nly a single unres o lved source 
at the position of S106 IRS4. ISolf Carse"^ (Il982l ) used 
the kinematics of the bipolar nebula to deduce that the in- 
clination angle of the system in the plane of the sky is ^ 75° . 

To determine the input parameters for our model of 
SI 06 we calculated a series of models with various Ms, Md, 
and Minfaii, where we fixed i, t's, ^'d, and cj. A good match to 
the Chandra count rate and the visual extinction, Av, was 
found for one particular model, with parameters as noted 
in Table [T] From this model we calculate a count rate (Ta- 
ble [3]) i n agreement with the 0.30±0.11 ks~^ Chandra de- 
tection (jGiardino et al.|[2004l V The A^H-star and A^ewc val- 
ues correspond to visual extinctions of — 22 and 39 
mag respectively, which are consistent with previously de- 



termined values (iGiardino et al.ll2004l ). The model param- 
eters for the stellar wind terminal velocity and mass-loss 
rate are consistent with an early B-type star. For the disk 
wind, the termina l velocity was set to the previous estimate 
of ^ 350 kms-^ (iDrew et al.lll993l ). whereas the mass-loss 
rate is roughly an order of magnitude lower than previously 
determined values of 1.6 - 1.8 x IQ-^Moyr" ^ (iFelli et alJ 
Il984l : iHoare et al.lll994l : iGibb k Hoa^l2007l ). This differ- 
ence may be due to wind clumping causing over-estimates 
in the mass- loss rate determinations from free- free emission. 
We note that amongst B-type stars, a mass-loss rate as high 
asl.8xlO~^M0yr~^ is consistent only with an evolved su- 
pergiant luminosity type, with main sequence stars having 
mass-loss rates typical ly two or more orders of magnitude 
lower re.g. |Priniai[l989h . 

4.2 Mon R2 IRS3 A 

The general structure of the observable nebula around 
IRS3 A is clearly monopolar, which suggests IRS3 A 
is embedded i n a d isk or torus with bipolar cavities. 
IPreibisch et al.l (|2002l l obtained high resolution (75 mas 
= 62 au) near-infrared H and K band images of Mon 
R2 IRS 3 which show a close triple system surrounded 
by strong diffuse nebulosity and three additional infrared 
sources within 3" of the brightest object IRS3 A {K ~ 7.9, 
M* = 12 — 15 Mq), which h as not yet developed an HII 
region (^Alvarez et al ] l2004bl ). For a review of the Mon 
R2 star-formin g regi on see ICarpenter k, HodappI (|2008h . 
IPreibisch et al.l (|2002h comment that the X-ray properties 
of IRS3 A and C cannot be explained by stellar wind mod- 
els where the X-rays are generated at shocks inherent to the 
winds as the emission is much h arder. This agrees with the 
comment by iKohno et al.] (|2002l ) that the observed plasma 
temperatures of ^ 4 keV (T > 5 x 10^ K) are ten times 
higher than plasma temperatures typical for early B-type 
main sequence stars. As has been shown in § (3] shock tem- 
peratures in a winds-cavity interaction can reach > 10^ K. 
In model Mon R2 (Table [T]) and the subsequent X- 



mined by I Alvarez et al.l (|2004al l from fits to observations 
with envelope radiative transfer models. The model count 
rate (Table [3)) agrees with the Chandra value of 0.166±0.041 
ks~^ determined by Preibisch et al. (2002). The stellar wind 
terminal velocity and mass- loss rate imply a late-O/early- 
B type central star. The mass infall rate for model Mon 
R2 is roug!:hly a facto r of 5 lower than that determined by 
I Alvarez et al.l (l2004al l. However, this difference is not sur- 
prising as the evolved state of the envelope in our model 
produces a dense layer of compressed cloud material along 
the cavity wall, which increases the column to the star. This 
feature is not present in envelope radiative transfer mod- 
els, which typically assume a smooth density distribution. 
As such, to attain a visual extinction which agrees with the 
observations we must decrease the mass-infall rate. 

4.3 AFGL 2591 



The - 16 M© (jvan der Tak k MentenI l2005l ) early B-type 
star AFGL 2591 sits at the centre of a powerful bipolar 
molec ular outflow cavity orientated in the e ast-west direc- 
tion (|Lada et al ] Il984l : iMitchell et aD Il992l ). The central 
source is fully obscured in the optical and J band. The exis- 
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Figure 14. Near-infrared speckle interferometry from 
iPreibisch et al.l (|2003) and contours from a Chandra X-ray 
image smoothed with a 1" gaussian super-imposed on the 
UKIDSS K band image of the AFGL 2591 region. The detected 
point-hke sources noted in Tableware labelled. 



tence of a large-scale torus perpendicular t o the outflow cav- 
ity is inferred from extended CS emission (|Yamashit a et all 
Nebulous loops are clearly seen in th e bisp ectrum 
speckle interferometry images of IPreibisch et alJ (|2003l ) . 
which are thought to mark the peripheries of o utflow bub- 
bles. W hen modelling high-resolution SO data, iBenz et al.l 
(|2007l l found the inclusion of X-ray emission (which may 
originate from outflow shocks) improved model flts. 

AFGL 2591 was observed by Chandra on 2006 February 
8 using the Advanced CCD Imaging Spectrometer (ACIS) 
S3 chip in very faint mode for an exposure time of 30.17 ks. 
Data were obtained from the Chandra Data Archiv^ and 
were processed in a standard way using t he Chandra In- 
terac tive Analysis of Observations fCIAO, [Fruscione et al.l 
l2006l ) software version 4.1.1. 

Point source detection was perfo rmed using the CIAO 
wavdetect task ([Freeman et al.l [20021 ) on the level 2 event 
flle. The process works by correlating the dataset with a 
number of "Mexican Hat" wavelet functions at different spa- 
tial scales to search for correlations. This process was per- 
formed at a significance threshold of 10~^, corresponding to 
less than 1 false source for the whole S3 field, and resulted 
in 4 sources detected in the vicinity of AFGL 2591 (Fig. [li] 
and Table 2]) and 62 sources in the whole S3 field. The point 
sources are identified with stellar point sources in a deep 
UKIDSS K band image of the region illustrated in Fig. [M] 

To correct for systematic errors in the Chandra aspect 
system (typically about 1 arcsec) , we searched for near- 
IR counterparts from t he UKIDSS (iLawrence et all [20071 ) 
Galactic Plane Survey ('Lucas et al.ll2008h catalogue across 
the entire CCD and identified bright counterparts for 12 of 
the detected X-ray point sources. These counterparts were 
used to correct the Chandra coordinates and register them 
to the UKIDSS reference frame. 

Source extra ction was performed using ACIS Extract 
(^AE jBroos et al.i>2002 ). an IDL-based package developed for 
ACIS data processing. The procedures used in AE are de- 

^ htt p : / / cxc . harvard . edu/cda/ 



Table 4. Detected point-like X-ray sources in the vicinity of 
AFGL2591. E is the median photon energy from the source. The 
net counts have been background subtracted. Pnot is the Poisson 
probability that the source is a chance coincidence of background 
photons — the "not-a-source" probability. 



No. 


RA 


Dec 


E 


Net 


Signif. 


Pnot 




(deg) 


(deg) 


(keV) 


counts 


(^) 




1 


307.3508 


40.1894 


3.9 


6.69 


1.77 


5.7e-8 


2 


307.3508 


40.1888 


3.2 


6.71 


1.77 


3.3e-8 


3 


307.3554 


40.1910 


2.9 


6.72 


1.78 


2.6e-8 


4 


307.3530 


40.1898 


4.3 


5.72 


1.59 


6.6e-7 



scribed in ? and ?. Model point spread functions (PSFs) 
were constructed for each source together with contours 
which contain >90% of source events. The background 
was estimated locally for each source using circular annuli 
around their respective PSFs that avoided the PSFs of other 
nearby sources. From the background subtracted net counts, 
the median energy, source significance and the Poisson "not- 
a-source" probability were calculated for each source. The 
results, source positions, and characteristics are listed in Ta- 
ble lU These sources are neighbouring low-mass stars and are 
not associated with the outfiow from AFGL 2591. 

AFGL 2591 was not detected as a point source, unlike 
S106 IRS4 and Mon R2 IRS3 A where the X-ray emission 
was unresolved. However, there is a tentative detection of 
diffuse emission (visible as the contours east of the central 
source in Fig. [T4)) with count rate of ^ 0.2 ks~^ from a cir- 
cular region, centered on AFGL 2591 with radius ^ 2.6^^ 
with a significance of 3.7 a ab ove the background level (de- 
termined using the method of IPease^e t al. 2006,, which ex- 
presses the probability that the observed counts cannot be 
explained by Poissonian variations in the background and 
is different to the signal-to-noise ratio of the extracted sig- 
nal). The origin of the observed emission is unclear and we 
interpret it as due to the interaction of the winds from the 
central object with the surrounding envelope. 

In our model of AFGL 2591 (Table [1]) we use a mass 
infall rate determined fr om a comparison of Eq. [3] with the 
Ivan der Tak et all (|l999l ) expression for number density and 
assuming the infalling cloud material to be molecular hy- 
drogen. We adopt a value of 15° for the cavity half opening 
angle fsee Ivan der Tak et allll999l : [Poelman van der TakI 
l2007l ). The stellar wind velocity is chosen to agree with 
an early B-type star. The disk wind velocity agrees well 
with the ~ 500 kms~^ velocities inferred from the ob- 
served high-velocity wi ngs of infrared absorption lines 
(I van der Tak et allfl999l ). Our model results show a count 
rate of 0.09 ks~^ which is consistent with the conservative 
upper- limit of 0.2 ks~^. The visual exctinction to the star 
(A y ^ 105 mag) agrees we ll with the 100 mag determined 
bv lvan der Tak et all (Il999l ) from JCMT observations of the 
C^^O J=2^1 and J=3^2 lines. The n iolecular gas with 
veloci ties of ^ 200 kms~^ observed bv Ivan der Tak et al.l 
(|l999l ) and interpreted as the entrainment of cloud material 
by the ionized wind outfiow can be readily explained by our 
winds-cavity model. 
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5 DISCUSSION 

For model Rl the efficiencies for converting tlie meclianical 
wind power into tliermal energy and tlien into X-ray emis- 
sion (= ivint/C) are :^ 5 x 10~^ and :^ 9 x 10~^ for tlie 
stellar and disk winds respectively, where ^ is the power of 
the specified wind. These values are representative of those 
calculated for all of the models with a common trend be- 
ing that the efficiency is roughly two orders of magnitude 
higher for the disk wind. In many cases the disk wind domi- 
nates the attenuated X-ray emission, despite the mechanical 
power being typically less than that from the faster stellar 
wind. However, the faster stellar wind is always responsible 
for the hardest X-ray emission. 

The model can provide a match to the observed X-ray 
characteristics of a sample of MYSOs. The intrinsic variabil- 
ity to the observed emission, caused by the dynamics and 
fluctuations in the position of the reverse shock, complicates 
placing constraints on model parameters. One possible av- 
enue for future work would be to attempt simultaneous fits 
to radio observations with the models. 

A recurring feature of the simulations was the delicate 
balance between the infall and the outflow. It proved im- 
possible to attain a set of parameters where the cavity wall 
remained stable and roughly stationary for t 2^ 10^ yr. The 
cause of this problem is the imbalance between the pressure 
in the infalling material and that in the outflows. Noting 
that the method used to describe the shape of the cavity 
wall is essentially arbitrary, a more realistic approach would 
be to have an initially fully embedded MYSO and to create 
a self- consistent cavity via a jet/outflow. 

Intuitively one would expect the opening angle of a cav- 
ity to be tied to the evolution of the central star(s), and that 
an initially narrow cavity resulting from a molecular outflow 
would be widened by outflows as the sta r evolves to finally be 



revealed as a main seque nce star (e.g . Velusamv k, Langerl 
1 19981 ). In a recent paper JCanto et all (| 20081 ) modelled out- 
flow cavities using a prescription with a time dependent 
opening-angle for the outflow. Although observational re- 
sults could be qualitatively reproduced, it is unclear from 
the wind-cavity models whether such an approach is physi- 
cally realistic. 

The model used in this paper considers the central ob- 
ject to be a single star. However, most massive stars form 
in binaries, which a llows the collision of winds between 
the stars as well (e.g. IStevens et al. l ll992l : IParkin k PittardI 
I2OO 8: Pittard 200^), and could p otentially generate consid- 
erably more X-ray emission (e.g. IParkin et al"]l2009l ). 



6 CONCLUSIONS 

Hydrodynamical simulations of the wind-cavity interaction 
around an embedded MYSO have been studied. The lati- 
tude dependent wind geometry we use incorporates both a 
stellar and disk wind. Using an extensive range of simula- 
tions we have examined the effect of varying different model 
parameters on the evolution of the cavity and the resulting 
observational characteristics (i.e. column densities. X-ray lu- 
minosities and count rates) . The main conclusions from this 
work are: 

• The collision of the winds against the cavity wall gener- 
ates a reverse shock (for cavity half opening angles uo < 60°) 



close to the star (< 500 au). The shock heated gas produces 
X-ray emission with an integrated count rate and spatial ex- 
tent in agreement with observations of MYSOs by Chandra. 
The position and shape of the reverse shock is dependent 
on the ram pressure in the inflow/outflow. Fluctuations in 
the position of the reverse shock cause variability of the ob- 
served emission on timescales of several hundred years, and 
possibly on shorter timescales which were not probed. 

• The amount of X-ray emission in the 4-10 keV band 
is dependent on the position of the reverse shock, which is 
strongly related to the stellar wind speed and the adopted 
mass inflow and outflow rates. 

• Integrated count rates in agreement with Chandra de- 
tections of MYSOs are obtained across a wide region of 
model parameter space, indicating that the generation of 
X-ray emission through the interaction of an outflow with 
infalling material is potentially a very robust process. 

• There is a limiting opening angle of the cavity for which 
a reverse shock is produced which resides close to the star 
and encloses the unshocked winds. For our adopted wind 
geometry we find this (full) opening angle to be 120°. 

• There appears to be a delicate hydrodynamic balance 
between the inflow and outflow. 

The model presented in this work provides a useful first 
insight into the interaction between winds from an MYSO 
and the surrounding envelope. As with many astrophysi- 
cal problems, there is a vast range of scales to consider 
in modelling the wind-cavity interaction, and the methods 
employed dictate the approximations that must be made. 
To accurately model the driving of the winds r equires high 
resolution in the wind acceleration region (e.g. IProga et al.l 
1998). Future models would benefit greatly from higher res- 
olution in the region close to the star and disk. This would 
allow the interaction between the outflow and inflow from 
the acceleration of the disk wind to the intersection point 
between the disk wind and the cavity to be examined. 
This would be an improvement on the current work as it 
would allow the stellar and disk winds to be self-consistently 
driven and the photoevaporation of the accret ion disk (e.g. 
iHollenbach et al.lll994l : [Richling k Yorkell2000h could be ex- 
amined. On larger scales, the evolution of the HII region 
around the star could be followed. 
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